import java.text.DecimalFormat;
import java.util.HashMap;
import java.util.Map;
import java.util.Random;
import java.text.*;
import java.util.*;
import java.io.*;

public class ThirdOrderMM extends MM {

  protected double[][][][] loge;

  private static double[][] toFirstOrder ( double[][][][] emat ) {
        double[][] aux = new double[emat.length][emat.length];
        for ( int i = 0 ; i < emat.length; i++ ) for ( int j = 0 ; j < emat.length ; j++ ) aux[i][j] = 0.0;
        for ( int i = 0 ; i < emat.length; i++ ) for ( int j = 0 ; j < emat.length ; j++ ) for ( int k = 0 ; k < emat.length ; k++ ) for ( int l = 0 ; l < emat.length ; l++ ) aux[k][l] += emat[i][j][k][l];
        return aux;
  }

  public ThirdOrderMM (Object[] esym, double[][][][] emat ) {
    super(esym,toFirstOrder(emat));
    for (int i=0; i<emat.length; i++) if (esym.length != emat[i].length) throw new IllegalArgumentException("MM: esym and emat disagree");
    this.esym = esym;
    nesym = esym.length;
    loge = new double[emat.length + 1][esym.length + 1][esym.length + 1][esym.length];
    for (int b=0; b < nesym; b++) {
      double fromstart = Math.log(1.0 / esym.length);
      loge[0][0][0][b] = fromstart;
      for (int j=0; j < emat.length; j++) for (int k=0; k < emat.length; k++) for (int l=0; l < emat.length; l++) {
	loge[j+1][0][0][b] = Double.NEGATIVE_INFINITY;
	loge[j+1][0][l+1][b] = Double.NEGATIVE_INFINITY;
	loge[j+1][k+1][0][b] = Double.NEGATIVE_INFINITY;
	loge[0][0][l+1][b] += emat[j][k][l][b];
	loge[0][k+1][l+1][b] += emat[j][k][l][b];
	loge[j+1][k+1][l+1][b] = Math.log(emat[j][k][l][b]);
      }
      for (int l=0; l < emat.length; l++) loge[0][0][l+1][b] = Math.log(loge[0][0][l+1][b]); 
      for (int k=0; k < emat.length; k++)  for (int l=0; l < emat.length; l++) loge[0][k+1][l+1][b] = Math.log(loge[0][k+1][l+1][b]); 
    }
  }

  public void print(PrintStream out) { 
    out.println("Symbol generation probabilities:");
    for (int b=0; b<nesym; b++) out.print(esym[b] + hdrpad);
    out.println();
    for (int i=1; i < loge.length; i++) for (int j=1; j < loge.length; j++) for (int k=1; k < loge.length; k++){
      for (int b=0; b<nesym; b++) out.print(fmtlog(loge[i][j][k][b]));
      out.println();
    }
  }

  public static MM train(Object[][] xs, Object[] esym ) {
    int nseqs  = xs.length;
    int nesym  = esym.length;
    double[][][][] E = new double[nesym][nesym][nesym][nesym];
    double[][][][] emat = new double[nesym][nesym][nesym][nesym];
    for (int i=0; i<nesym ; i++) for (int j=0; j<nesym ; j++) for (int k=0; k<nesym ; k++) for (int l=0; l<nesym ; l++) E[i][j][k][l] = 0.0;
    ThirdOrderMM mm = new ThirdOrderMM(esym,E);
    for (int s=0; s < nseqs; s++) {
        Object[] x = xs[s];
        for (int i=0; i < x.length; i++) for (int j=1; j < x.length; j++) for (int k=2; k < x.length ; k++) for (int l=3; l < x.length ; l++) E[mm.getSym(x[i])][mm.getSym(x[j])][mm.getSym(x[k])][mm.getSym(x[l])] += 1.0;
    }
    double count = 0;
    double countZero = 0;
    for (int j=0; j<nesym; j++) {
        double Eksum = 0;
        for (int k=0; k<nesym; k++) for (int l=0; l<nesym; l++) for (int b=0; b<nesym; b++) {
		count++;
		if ( E[j][k][l][b] == 0.0 ) countZero++;
		Eksum += E[j][k][l][b];
	}
        for (int k=0; k<nesym; k++) for (int l=0; l<nesym; l++) for (int b=0; b<nesym; b++) emat[j][k][l][b] = ( E[j][k][l][b] + 1.0 ) / ( Eksum + 0.0 + ( nesym * nesym * nesym ) );
    }
    if ( debugTrain ) System.err.println("Model sparsness : " + countZero + " in " + count + " ( " + ( countZero / count ) + " % )");
    return new ThirdOrderMM(esym,emat);
  }

}
